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Abstract 

This report covers an analytical and computer 
study of the dynamic energy equations that describe 
the physical phenomena that occurs in a Stirling cycle 
engine. The basic problem is set up in terms of a set 
of hyperbolic partial differential equations. The 
characteristic lines are determined. The equations 
are then transformed to ordinary differential equations 
that are valid along characteristic lines. Computer 
programs to solve the differential equations and to 
plot pertinent factors are described. 
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INTRODUCTION 

The dynamic energy equations treated in this report describe the 
energy conversion phenomena that occurs in a Stirling cycle engine. 

Analytical and numerical techniques that are somewhat different from those 
used by previous investigators have been applied. The problem is rather 
complex in that it involves the non-steady oscillatory flow with compres- 
sion and expansion of a gas in a heater, regenerator and cooler configuration. 

The basic problem for this study has been an analytical and numerical 
description of the energy conversion in a Stirling cycle engine by means 
of Dynamic Energy Equations, and a definition of the temperature, density 
£uid pressure as these variables change with time and with position during 
the cyclic operation of the machine. 


PREVIOUS WORK IN THE FIELD 

There is an extensive body of literature in the field of Stirling cycle 

machine analysis and design. Only a few of the basic references are listed 

1 2 
in this paper. Walker is a standard basic reference. Martini provides a 

3 4 5 

concise Stirling cycle engine design manual. Finkelsteln , Uriell , Martini , 

6 7 8 

Tew, Jeffries and Miao , Tew , and Martini present various analytical and 

computer approaches. Rios^, Urieli^, and Kirkley** are Ph.D. thesis on 

12 13 

various aspects of the Stirling cycle. Larson and Larson cover some of 

the analytical concepts and computation techniques discussed in this report. 

14 

Lorenzo and Daniele presents aspects of Stirling engine performance and 
control. The above references include ideal Schmidt cycle analysis, lumped 
control volume analysis, regenerator effectiveness studies, nodal analysis, 
and general thermodynamic studies. The authors cited have done a great deal 
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of Stirling cycle work and have written many other papers. Sons of the 
cited references contain extensive bibliographies. The other references 
listed present mathematical, physical and computational information that 
is pertinent to this study. These latter references are cited at appro- 
priate points in following pages of this report. 


THE PHYSICAL PHENOMENA 

The working gas in a Stirling cycle engine is sealed between a com- 
pression piston and an expansion piston. The gas oscillates through a 
heater , a regenerator and a cooler as the pistons reciprocate. In Stirling 
machines the expension piston reciprocates with t phase difference from 
the compreselon piston so the gas is periodically compressed and expanded. 
Heating and cooling ia provided so that power is produced from heat in an 
engine configuration (or heat energy produced from work for the heat pump 
configuration). The gas flow region for s typical configuration la illus- 
trated schematically in Figure 1. The regenerator is a passive unit 
cyclically absorbing heat from, and releasing heat to, the working gas. 

The gas flow ia driven by the expansion and compression pistons which move 
with an approximate sinusoidal motion with a phase difference (usually about 
90*). In this study sinusoidal motion with a phase difference of 90 degrees 
has been assumed, with temperatures and phase relationship typical for an 
engine. A change in input, temperatures, and piston phase will provide heat 
pump analysis. 

The cyclic phenomena occurs at the frequency of engine rotation. 

Under the quasi steady state conditions the gas conditions will repeat 


9 

f 

every 360° of piston motion. The gas boundary is assumed to be at the 
respective constant wall temperatures in the heater and cooler. In the 
regnerator the mesh temperature varies with time and position. The time 
variation is cyclic. The position variation covers the range from near 
the cooler temperature at one regenerator boundary to near the heater 
temperature at the other end. 

ENGINEERING APPROACH TO THE PROBLEM 

The physical phenomena that occurs in a Stirling cycle engine can be 
described by a set of dynamic energy equations. In a following section of 
this report it is shown that those Dynamic Energy Equations, a set of partial 
differential equations, are of the hyperbolic type. The specification and 
numerical solution of the descriptive equations is the object of this study. 

The engineering approach to the problem includes three major aspects: 

1. A mathematical description of the phenomena that occurs with the 
oscillation of the gas. This description includes appropriate formu- 
lation of the laws for conservation of mass, conservation of energy, 
conservation of momentum and the equation of state for gas. Convective 
energy transfer is a major factor in the problem. 

2. Formulation of the differential equations and analytical expressions 
in a form appropriate for digital computer solution. An Important 
aspect of the engineering approach has been to formulate expressions 
that describe applicable conditions with sufficient accuracy for the 
problem and that will result in a tractable overall problem statement. 
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3. Development and eelectlon of suitable computation techniques. 

A fourth-fifth order Runge-Kutta numerical solution technique has 
been adapted to solve the ordinary differential equations that result 
when the characteristic directions of the hyperbolic system have been 
determined. The computation technique Includes a subroutine to compute 
appropriate segment lengths for the cooler, regenerator and heater. 

The analytical, numerical and logical techniques are discussed in 
following sections of this report. 

INNOVATIONS IN THIS STUDY 

This study Includes a number of techniques which, to the author's 
knowledge, have not been applied previously to the analysis of Stirling 
cycle machines. These include the partial differential equation classi- 
fication with a characteristic equation approach to the problem; treatment 
of cooler, regenerator and heater as connected units with selectable 
transition regions (instead of a more conventional lumped parameter approach) ; 
and the use of heater and cooler increments that grow at a constant ratio as 
distance from the regenerator increases. 

The Runge-Kutta integration program (RX745) was selected from the 
literature^ and used with minor modifications for this study. The gas 
velocity is defined in terms of a linear Lagrange interpolation to define 
its velocity between the compression and expansion pistons. This provides 
an analytical expression for the gas velocity as well as analytical expres- 
sions for the velocity derivatives with respect to poslMon and with respect 
to time. Details of these techniques are presented in following sections of 
this report. 
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The simultaneous partial differential equations are shown to be of the 
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hyperbolic type. Since the gas velocity is small relative to the sonic 
velocity in the gas, the wave phenomena can be neglected. It should be 
noted that the equations are of the hyperbolic type when wave phenomena is 
neglected. Thus the characteristics of this analysis are not the charac- 
teristic Mach lines of compressible flow. The characteristic lines, of 
variable slope, are functions of position and of time. The slope of the 
characteristics is dependent on the direction and velocity of the gas flow 
in the system components. The analytical expression for the gas velocity 
provides a good approximation to the actual gas velocity and its derivatives 
in the heater, regenerator and cooler and provides the necessary Information 
to define the varying characteristic directions. 

The computation takes advantage of the fact that along the character- 
istic lines the position variable is eliminated, thus transforming the partial 
differential equations into a set of simultaneous ordinary differential 
equations. The resulting ordinary differential equations are solved numeri- 
cally using a fourth-fifth order Runge-Kutta integration technique. The 
computer program includes the determination of the characteristic directions 
that apply for the current step and an automatic computer adjustment of the 
step size. The Runge-Kutta Integration along the characteristic lines 
provides a solution technique that is more stable than finite difference 
techniques frequently utilized in the solution of partial differential 
equations. An initial cosine curve assumption for the regenerator mesh 
temperature quickly produces a cyclic repetition of the computed gap. 
temperature, density and pressure conditions. The temperature, the variable 
of primary Interest, is automatically plotted using (off-line) computer 
graphics equipment. 


- 5 - 


THE BASIC EQUATIONS 

The energy transfer in s Stirling cycle machine can be described by 
three basic principles: 1) conservation of mass , 2) conservation of energy, 

and 3) conservation of moment us. The general equations on a unit volume 
basis ara 

3p /3t + 7 .(p O ) - 0 (1) 

© o o 

3/at(p g C v T g )/at - - V • [p g U g (h c +U g U g / 2GJ)+K g (VT g ) ] + aw/at - “ c A c (T g -T B ) (2) 
pgl (u g • vu g ) + au g /atj + 7? g - F y (3) 

The momentum equation (3) will not be explicitly required In the set 
of simultaneous equations if an analytical expression for the gas velocity 
is available. In a following section of the paper an analytical approxi- 
mation for the gas velocity is derived. It should be noted that while 
equation (3) is not utilized, the momentum factor is Implicitly carried 
along in the computation through the gas density, velocity and pressure 
terms. 

The enthalpy of the gas is typically much greater than the kinetic 
energy of the gas. Also the wave velocity in the gas is normally much 
greater than the gas velocity; therefore, the times of Interest for the 
oscillatory flow are such greater than pressure wave propagation times. 
Assuming these fluid flow characteristics, the kinetic energy terms and the 
second derivative terms in equation (2) can be neglected. It should be 
noted that the compressibility of the gas is taken into account when . 
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ORIGINAL FAukJ id 
OF POOR QUALITY 

computing volume, pressure, density and temperature levels and changes in 
these parameters. 

In all following equations it is assumed that the flow is one-dimensional. 
Therefore the vector notation will be dropped. (For one-dimensional flow the 
differential operator 7 is a scalar operator, 7 - 3/3x). 

With the above assumptions, the partial differential equations defining 
the temperature and the density of the gas can be written in the matrix form 


—■ — 





U 1 0 0 


p 

. 

- p u 

g 


X 


g g* 

r Y« T , wB . p . 


p t 


- yp T U - h A (T -T )/p C 
,H g g gx ccgm g v 



T 

X 


mm. 


T 



t 



The subscripts x and t on the p, T and U terms indicate partial 
differentiation with respect to x and t respectively. Thus, for example, 

0 t "3p/3t and U^-SU^/ax etc. 

In addition to the equations describing gas conditions, an equation 
describing the regenerator material temperature is required. 

0» C »< 3 V 3t > * <T -T ) + K 7 2 T (5) 

mm m ccgm mm 

In this equation the subscript m refers to the heat transfer material. The 
temperature of the material receiving and providing heat in the cooler and 
in the heater will be assumed to be constant. 
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AN ANALYTICAL APPROXIMATION FOR THE CAS T 

VELOCITY AND ITS DERIVATIVE WITH RESPECT TO X 

A Stirling cycle engine will be teken ee the specific type of aechine 

in this peper. S lapis changes in nomenclature and detail would permit 

generalisation to othev Stirling aechines. The portion of e Stirling 

cycle engine of interest is illustrated schematically in Figure 1. 


Compression 

Cylinder 



Expansion 

Cylinder 




j»Cooler -*j 


f 




U- Regen J 


Heater 



■4low Transition 


(Dead Volume) 

Figure 1. Schematic Diagram of Stirling Cycle Engine 


The velocity of the gas is a function of the velocity of the compres- 
sion and expansion pistons , the flow geometry* and the flow resistance. The 
pistons will be assumed to nova with simple sinosoidal motion with the 
expansion piston lagging the compression piston by 90* 


U ■ R wsinwt 
c c 


( 6 ) 


U - R oislnftdt - 90*) (7) 

e c 

An "ideal” gas velocity would occur with sero flow resistance. This 
is equivalent to assuming uniform gas pressure in the region between the 
pistons (flow pressure drop and pressure waves are neglected) * and also 
assuming one-dimensional flow (even though the flow area changes). The gas 
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velocity between the pistons cen be detetermined by s linear Lagrange 
interpolation. Note that if the flow area were uniform and equal to the 
piston area the ideal gas velocity would be 


U . " < u „ x . *■ u „0 / < x „ ♦ x > <®> 

gi c e a c c e 

where x end x are the distances from the point of Interest to the expan- 
c e 

sion and compression pistons respectively. The geometry of a typical flow 
path between the pistons is rather complex. In that case the distances x £ 
and x # from the pistons are not as significant as the gas volumes in deter- 
mining the ideal gas velocity. 

With the varying flow area generalizing equation (8) the ideal gas 
velocity (a function of x and t) will be approximately given by a "volume" 
Lagrange interpolation. 

U . " <V A f > < u r v . + / < v r ♦ V < 9 > 

where and are the actual gas volumes from the point of interest, x, 
to the compression and to the expansion pistons respectively. 

The ideal velocity flow corresponds to nonvlscous or zero pressure 
drop flow. The actual flow results in a pressure drop. That pressure drop 
depends on the length of the flow path, the fluid velocity, the flow path 
configuration, and the Reynolds number. The periodic flow will result in 
a flow pressure drop that is a function of time and position. 

For the present application it will be convenient to define the effect 
of flow resistance on the gas velocity. So that 


U 


8 



( 10 ) 
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Consider a point x along the fl'v path. To be specific assume x la 
In the regenerator. Because of flow resistance the pressure and the density 
will be greeter at points In the direction of the piston "driving" the flow 
than It will be at points In the direction of the receding piston. The flow 
reslstance-pressure-danslty characteristic will result in a gas velocity that 
Is less In absolute value than the Ideal at all points alcr ; the flow path of 
the oscillating fluid expect at the ends In contact with a piston. 

The pressure drop can be expressed 

BP - fp g LU g */2CD (11) 

The value of the friction factor (f) for laminar flow Is 


f - 64/Re - 64u/( Pi |0 g |0) (12) 

7or turbulent flow the friction factor is dependent on the Reynolds nuaber 
(Re) and on a tube roughness factor. Typically Re will be less than 15000. 

For Re values between 2100 cod 15000 the friction factor can be approximated 
by the equation 

t - C lf ♦ C 2£ /Re (13) 

C 1£ and C 2£ are taken to approximate a selected curve. Kota that In the 
laminar regime (Re < 2100) C 1£ - 0.0 and C^ f - 64.0 

18 

From values tabulated In reference It appears that the viscosity (y) for 
hydrogen and for helium can be ipproxlaated by linear functions of 
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temperature (for temperatures in the range 300 < T < 1100°K). With this 

© 

approximation 

U ■ a + bT (14) 

S 

Using equation (13) in equation (11) and taking acceleration effects into 

account by adding the 3U /3t and 3U /3x terms* gives 

© 6 

AP - - K p (C lf 4C 2f /Re)p* f |U g |U g /(2GD f ) - p(U gt +D g U gx ) x f /G (15) 

Considering the fluid as a spring, the force AP A, to produce the flow 

© 

is obtained by a "deflection" or decrease in x, say Ax. Let K be a "spring 

s 

constant" then* for small Ax 

£P A. • K Ax (16) 

g t s 

The change in pressure from the uniform pressure (zero viscosity) case 
is given by a pair of integrals, one with limits from the compression piston 
to x (the point of interest) * and the other with limits from x to the 
expansion piston. That pressure effect can then be interpreted in terms of 
its effect on the position in equation (16). 

It is reasonable to neglect any change in temperature due to the rela- 
tively small flow pressure drop because the cooler* regenerator and heater 
are designed for high heat transfer. Then noting the total pressure drop 
along the flow path has a contribution from the point x to the compression 
piston* and one from x to the expansion piston 

|AP C /P| + |AP e /P| - lAV^/Vj +|AV e /vJ (17) 
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Then letting A^Ax » AV fi • AV # , and AP “ AP fi + AP^ 


l«l •!»( M <V v . )/v c v . 


a»> 


(The absolute values are used because the effect of the sign of AP and 
Ax will be taken into account in equations (20) and (21).) 

The change in distance (Ax) in equations (16) and (18) results in a 
change in gas velocity (AU^). As suae that 

AO/U • KAx/x. (19) 

g g t 

Then equating the pressure drop as given by equations (IS) and (18), 
inserting Ax from (19) and solving for AD (and designating the constant 

o 

term as K y ) gives 


- K v (C lf 4C 2f /R.)(, g X f V c T e U 1 gl /(2GD { P g (V c+ V t ) J ) (20) 

The partial derivative of U with respect to x is approxinately 

*0, A p <VV 11 • VcVf (3C lf p0 *« + 2C 2f‘ , l D »il >] + (21) 

3x “ (V c+ V.) 2GD f P g (V c +V e ) iI 0 { 

- VVVV^ tC H +C 2f / ** 1 D ’*l 

2eD t T g (v c +vJ 2 
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Typical velocity curves are illustrated by Figure 2. The effect 
of the flow resistance on the gas velocity is exaggerated by the curves 
on the figure. The effect is most apparent in the slightly curved 
lines for 150° and for 330°. 



FIGURE 2 TYPICAL GAS VELOCITY PATTERN 






CLASSIFICATION OF THE 
PARTIAL DIFFERENTIAL EQUATIONS 
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The analytical approximations to U and 3U / 3x will be used in the 

© © 

sequel to simplify the treatment of the dynamic energy equations. The gas 
density and temperature are functions of position and of time. The total 
differentials of these terms are 


dp - (3p_/3x)dx + (3p /3t)dt 

o © o 

(22) 

dT - (3T /3x)dx + (3T /3t)dt 
8 8 8 

(23) 


Appending these equations to equation (4) 


U 

8 

1 

0 

T Vg 

T 8 

YP ® 

T 8 8 

d* 

dt 

0 

0 

0 

dx 







0 


p x 

- 

“ P 8 U 8* 

P 8 


p t 


- VsVWVV'V* 

0 


T 

X 



dt 


T t 


dT 

8 


(24) 


This system of partial differential equations is of the hyperbolic* 
parabolic* or elliptic type when the discriminant of the characteristic 
equation is positive* zero* or negative respectively. (Ames reference 
The characteristic equation of this set of simultaneous partial 
differential equations is obtained by setting the determinant of the 
coefficient matrix equal to zero. The characteristic equation is 


- YPgUg 2 dt 2 + p g (y+l)U g dtdx - Pgdx 2 ■ 0 


(25) 
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Dividing through by p dx 2 (p is nonzero) gives a quadratic equation in 

s s 

dx/dt. The discriminant of that equation is 

B 2 - 4AC - U a (y-l) 2 (26) 

© 

Note that this discriminant is greater than zero (except locally in time 
and space when U is zero. Therefore the set of equations is hyperbolic 

O 

(except parabolic when U ■ 0) . 

8 

CHARACTERISTIC DIRECTIONS 

Solving the characteristic equation for dt/dx provides the charac- 
teristic directions. Then assigning ot and $ to the distinct directions 


dt/dx) - 1/U 
o g 

(27) 

dt/dx) 6 - l/ Y U g 

(28) 


Figure 3 is a computer drawn plot of a characteristic line network 
in the regenerator portion of the flow path. 



FIGURE 3 PLOT OF ALPHA AND BETA CHARACTERISTIC LINES 
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DERIVATION OF DIFFERENTIAL EQUATIONS OF POOR QUALITY 

ALONG THE CHARACTERISTICS 

Along the characteristic lines the determinant la aero. The equations 
will have a solution only if the determinant is also sero when the right 
hand column vector la substituted for any column of the matrix. Then 
setting up the equation in terms of dt/dx and successively inserting the 
a and 0 characteristic values for dt/dx will reduce the partial differential 
equation to a set of ordinary differential equations. It is necessary to 
select appropriate columns for replacement by the right hand vector ao as 
to get an ordinary differential equation for each of the dependent variables. 
(In some cases a trivial result may occur.) 

Substituting the right hand column vector for the first column gives a 
determinant that must equal zero for a solution. 


- p U 
8 8 * 


Yp T U + 

8 g gx 

- h A (T -T )/C p 
c c g m' v g 


W ."« 


dp g dt 0 

dT 0 dx 

8 




(29) 


After substitution of dt/dx) ■ 1/U this determinant reduces to the 

a g 

equation 

(p g 3U g /3x + dpg/dt) ( Y p g U g - p g U g ) - 0 (30) 


Since y f 1 1 the first term must equal zero. Therefore the density is 
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defined along the o characteristic by the ordinary differential equation 


dp/dt - - p(3U/ 9x) 

8 ga 8 


Along the 6 characteristic the substitution dt/dx)^ ■ 1/yUg gives 


(p g 3Ug/3x + dp g /dt) (YPgU g - YP g U g ) ■ 0 


The second term is identically zero so another substitution is 
required. Substitution of the right column for the second column gives 
the same result as above* i.e., defines dp/dt along the a characteristic* 

O 

but does not define dT /dt. 

8 

Replacing either the third or fourth column of the matrix by the right 
hand vector gives an equation including dT /dt and dp/dt. Taking the 

O © 

value for dp /dt along the a characteristic 
8 


dp/dt - - p (3U /3x) 

g a g 


and using the $ value* dt/dx) - 1/yU provides an ordinary differential 

p 8 

equation defining the gas temperature 


dT /dt - - T. a U a (v - p U /p a U J - h A (T a -T )/C p a 
g gS gxB 1 a gxa 8 gxg c c gB nr v 8 


In this equation the a and $ subscripts Indicate that the various functions 
are to be evaluated along the a and 8 characteristics respectively* thus 

U gxa 3 V 3X> a 


( 31 ) 


( 32 ) 


( 33 ) 


( 34 ) 


( 35 ) 
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THE COMPUTER PROGRAM 


The computer program CD EE (characteristic dynamic energy equations) 
ia designed to aolve equations (3), (5), (33) and (34) using equations 
(27) and (28) to provide appropriate characteristic points* using 
equations (9)» (10)* (20), and (21) to provide gas velocity information. 

The gas velocity functions uncouple the momentum equation (3) from the 
system of equations. Mote that equations (33) and (34) do not explicitly 
involve position* and include position derivatives only in the dU^/dx terms. 
This latter term la provided by the analytical expression (21). Many other 
mathematical equations and logical criteria are also utilized in the program. 
The Runge-Kutta RKF45 integration program utilizes computer-determined 
basic time increments (H) with weighted integration points at H/4* 3H/8* 

H/2, 12H/13 and H. It is accurate to fifth order (sixth order locally) 
according to reference* 7 . The characteristic directions are utilized to 
determine appropriate "characteristic points" for evaluation of the 
variables. 

The utilization of the characteristic directions and the analytical 
expression for tha gaa velocity provide a significant advantage for numerical 
solution of the dynamic energy equations. The advantage is that the usual 
finite difference solution of the partial differential equations can be 
replaced by a rather sophisticated Runge-Kutta integration routine, des- 
cribed in detail in reference* 7 . Adams type predictor correction integration 
techniques could be used but have not been applied. 

The major parts of the program are: 
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A. CDEE, Characteristic Dynamic Energy Equations 

1. The MAIN program 

2. The CHIC subroutine 

3. The REGEN subroutine 

4. The Runge-Kutta Integration program RKF45 
with subroutines* RKFS, and FEHL 

B. DEEPP, Dynamic Energy Equation Plotting Program 
The MAIN Program of CDEE 

The main program provides general program control. This portion of 
the program provides for specification of the computation parameters* and 
logical and physical data. MAIN provides for specification and computation 
of derived initial and boundary conditions. The variables used in the 
program are defined by "consent" statements and are provided with default 
values. MAIN provides for all input and output. Input variations are entered 
through three NAMELISTs so that no special input FORMAT statements are required. 
In practice the NAMELISTs are used to change the default values to those 
desired for a specific analysis. 

The program output is from MAIN. It contains all WRITE statements and 
all FORMAT statements. Error messages are handled by setting flags so that 
such output is also from MAIN. 

The variables and parameters used in the program can be divided into 
the following four categories: 

1. Computational data that define time and position increments* 
accuracy requirements, the extent of the computation, and 
other program control and bookkeeping functions 

2. Logical data that can be set up to permit alternate conditions* 
alternate fluids and differences in operational or computational 
logic 
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3. Numerical data for the independent variables and parameters that 
define the physical configuration, the material properties and 
the operating conditions for a particular analysis 

4. Dependent variables or computed values that comprise the computer 
output and provide the desired Information on system characteristics 

The first three groups of variables and parameters are included in three 
NAMELlSTs: COMDAT, LOGDAT, and NUMDAT, respectively. Default values for all 

inputs are provided in MAIN so that only those requiring different values 
need be entered through the appropriate NAMELIST. 

The computational data is listed alphabetically In a NAMELIST called 
COMDAT. This data Includes such Items as number of Increments to be used in 
each part of the computation. Initial and final time, or degrees of rotation. 
(The actual time increment used varies during the computation and is auto- 
matically adjusted in RKFS, the Runge-Kutta integration subroutine.) 

The program provides for selection of different working gases, alter- 
native computational logic, and various assumptions through the specification 
of selected logical data which are listed alf habetically in a NAMELIST called 
LOGDAT. 

The numerical data and Independent variables for the program are 
listed in a NAMELIST called NUMDAT. This includes data that describes the 
physical configuration, various material properties, initial conditions, 
boundary conditions, and the physical constants that are required for the 
computation. 

All NAMELlSTs are arranged alphabetically to facilitate checking. 
Virtually all Input data is Included in the NAMELlSTs. However, because of 
the NAMELIST input characteristic (and because default values are provided) , 
only specific Items In the list to be changed need be typed into an input 



NAMELIST. All the NAMELISTs are read by MAIN before execution and all are 
printed out In the standard free format at the beginning of the program 
execution, and also at one point aljng with the formated output. This 
provides a convenient check of the actual parameter values that were used 
for a particular computation. 

The MAIN program Includes segments to insure that total mass is con- 
served, to compute work output, and to enter data on tape for off-line 
plotting by DEEFP. 

The CHIC Subroutine 

The regenerator la the region of greatest Interest. It has by far the 
greatest temperature variation, ranging from the temperature of the cooler 
to the temperature of the heater. The uniform Increment length In the 
regenerator (LR/NINCR) Is computed from the selected number of Increments. 
Variable Increment lengths are uaed In the heater and cooler to permit a 
reasonable total number of increments and still include all of the long, 
relatively uniform temperature heater and cooler. The subroutine CHIC 
(Cooler, Heater Increment Computation) is provided to adjust the size of the 
Increments in a progressive manner. The increment on either side of the 
regenerator is the same length as those In the regenerator. A ratio (R) Is 
computed such that a selected number of Increments will cover the full length 
at the heater (or cooler). The length of the Ith Increment (in the heater, 
for example) is 

XINC(I) - XINCR * R**(I-1) (36) 

The ratio R is computed from the equation 

XINCR (1 + R + R 2 +...+ R N-1 ) - LH (37) 
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The number of Increments In the cooler, regenerator, end heater are 
specified independently in the computation data NAMELIST (COMDAT) . 

The printout of pertinent data in a controlled format at specified 
time (or angular rotation) increments and at specified evaluation positions 
is provided by MAIN. 

The REGEN Subroutine 

REGEN is the portion of CL EE In which the gas velocity is computed and 
the equations to be Integrated are evaluated. This subroutine provides for 
computation of required characteristic points which do not in general 
coincide with the specific position increments into which the system Is 
divided. REGEN utilizes a number of logical variables that permit a variety 
of computational assumptions. They include use of hydrogen or helium, 
various fluid flow area change configurations, and other flow character* 
lstlcs. The area variations include linear area variation (LAV), linear 
velocity variation (LW) , and uniform velocity ration (UVR) at successive 
evaluation points. The primary function of REGEN is evaluation of the gas 
velocity and evaluation of the time derivatives of the density, the gas 
temperature, and regenerator mesh temperature. REGEN is used as an EXTERNAL 
function which is called by arguments in the RKF45, RKFS and FEHL subroutines. 

The RKF45 Runge-Kutta Integration Program 

The Runge-Kutta integration is carried out by three subroutines 
described below. These subroutines, described in detail in reference^ were 
modified to permit computation using the varying network of characteristic 
points in place of the usual regular mesh network. 
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RKF45 is an Interface subroutine that implements the fourth-fifth 
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order Runge-Kutta integration of a system of first order ordinary 
differential equations. This subroutine is the bookkeeping portion of the 
overall Integration program which Includes RKFS and FEHL in addition to 
RKF45. 

RKFS is a control subroutine that checks accuracy requirements and 
provides for automatic adjustment of the time step. The actual numerical 
integration is carried out by the subroutine FEHL which is called from 
RKFS. 

FEHL Integrates the system of first order equations (in REGEN) . Initial 
values and the Initial derivatives are specified at the starting time. 
(Initial in this context means simply the values at the start of a particular 
step.) FEHL advances the solution over one time step that is set by RKFS. 

The variable time step characteristic of the program is provided by RKFS so 
that FEHL receives a time step that is fixed (but may be different at 
another call) . 

DEEPP, a Plotting Program 

The gas temperatures computed by CDEE are written into a library file. 

A particular library file called CEDATnnn is retrieved for off-line plotting 
by computer graphics equipment. The plotting program DEEPP (Dynamic Energy 
Equation Plotting Program) is an interactive program written in BASIC. The 
program is flexible so that it can be run with temperatures in degrees 
Rankine or in degrees Kelvin; the scale can be adjusted for the size plot 
desired. The plot can be either temperature versus time at a specified 
position or temperature versus position with crank rotation as a parameter. 
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Typically the cooler and heater are ouch longer than the regenerator* but 
the temperature variation In the regenerator la of greater lntereat than 
the relatively uniform cooler and heater temperature. To provide a good 
"picture" of the gaa temperature profiles In the regenerator the major 
portion of the plot (actually the center half) la devoted to regenerator 
gaa temperature. The program la aet up to uae different scales In the 
three regions* or fer a true scale but with only a portion of the heater 
and cooler Included. In the latter case the scale is set by the regner- 
ator, which Is half of the plotted abscissa* the part of the heater and 
cooler Included on the plot is equal to half of the regenerator length. 

This BASIC program Is Interactive. The various alternatives are obtained 
by operator response to questions printed on the screen. 

Computational Parameters 

The numerical values associated with a particular physical configuration 
are very significant In an analysis. The determination of appropriate 
values Is usually a routine computation, however occasionally the values 
and the effects are subject to interpretation. In some cases significant 
research efforts may be required to determine the effects. 

The ratio of heat transfer area tc flow volume (or perimeter to flow 
cross section) Is a very significant parameter for the cooler* whs heater 
and the regenerator. Por units with simple round tubes the ratio Is simply 
4/dlam. with an appropriate dimension. For a regenerator comprised of many 
layers of fine stainless steel screen the effective value is not obvious. 

For this study the following values were used: 
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A * 3700 m 

cc 

A - 32000 ra 
cr 

A ■ 1300 m 

cn 


The heat transfer was defined In terms of the Nusselt number (Nu - hd/k) . 
Laminar flow was assumed when the Reynolds number (Re * pUd/y) Is less 
than 2100. It Is assumed that the flow is turbulent when Re Is greater 
than 2100. The following equations were utilized for convective heat 
transfer: 

Nu - 3.65 for laminar flow 

Nu - 0.023 (Re)0‘®(Pr)0.4 for turbulent flow 


In the latter equation the Prandtl number (Pr * C^y/k) was taken as 

constant (Pr « 0.70 for hydrogen and Pr ■ 0.73 for helium). The gas 

conductivity k and the gas vlscc^ty (y) were approximated by linear 

functions of gas temperature. These approximations provided values reasonably 

18 

dose to values as given in the handbook for the limited temperature range 
that is applicable for this study. 

The pressure drop effect on gas velocity is assumed to vary with 
distance between the compression and the expansion pistons according to 
the relationship as given by equation (9) , repeated here for convenience 



(A/A f )(U V + U V ) / (V + V ) 
pz ce ec c e 


Note that at the compression and at the expansion pistons the gas velocity 
will be equal to the respective piston velocity as it must be. Between 
the pistons the gas velocity is a gas volume Lagrange average. 
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Using the equation (9) and the assumptions pointed out previously, 
the change in gas velocity from the Ideal Is given by equation (20) . Note 
that this expression (that Is the assumptions leading to this expression) 
results In a gas volume effect that takes the gas volume at each end Into 
account. The maximum effect of this factor Is at the center as would be 
expected. That Is, the maximum decrease In gas velocity occurs at the mid 
point between the pistons (In terms of volume) . The maximum value of the 
V c V e /(V c +V ft ) 2 term is 1/4. To bring the net gas velocity effect of the flow 
resistance up to the assumed value .a pressure increment coefficient (PICO) 

Is introduced as a multiplicative factor In the pressure drop and in the 
velocity drop equations. One additional factor Is Introduced because it 
appears that the pressure drop, particularly that in the regenerator. Is 
understated by the classical pressure drop expression (equation 11). This 
seemed to be particularly significant for the helium cases. To compensate 
for this perceived effect a factor that will provide an appropriate nonlinear 
gas velocity effect (with greatest significance at the higher gas velocities) 
is needed. This has been provided by an Increase In an exponent (PIEX) on 
the absolute gas velocity from 1.0 to PIEX with a value greater than one, 
typically a value of 1.25 has been used. PICO and PIEX are Included in the 
NUMDAT NAMELIST so values other than the default values can be used when that 
is deemed advisable. 

RESULTS 

The variation of gas temperature with position and with time Is one of 
the more distinctive characteristics of Stirling cycle engines. Typically 
temperature varies within a range of ± 30° C of the cooler temperature In 
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the regions close to the compression cylinder. In the heater typical 
variations are within ± 130® C of the heater temperature. (Specific 
typical values are given below.) This variation is due to the oscillation 
of the gas between the various components along the flow path as the gas is 
subject to periodic compression and expansion. It should be noted that 
typical volumes of the cooler, regenerator, heater and dead spaces are 
such that the gas does not traverse the complete path. For example. In 
the GPU-3 engine the total volume between the cylinder is approximately 
1.6 times the total displacement volume of one cylinder. The relative 
volumes (with values comparable to those for a GPU-3) based on the total 
maximum volume, which for the assumed configuration occurs when the com- 
pression piston crank is at 135°, are as follows: 

At 135° At 315° At At 



(Max Vol) 

(Min Vol) 

TDC BDC 

Compression cylinder 

25.51 

4.4%* 

0%* 

30%* 

Dead volume cooler region 

11.4 

- Same as 

at left 

- 

Cooler 

2.6 

- Same as 

at left 

- 

Regenerator 

5.3 

- Same as 

at left 

- 

Heater 

11.3 

- Same as 

at left 

- 

Dead volume heater region 

18.4 

- Same as 

at left 

- 

Expansion cylinder 

25.5 

4.4 * 

15 * 

15 * 

Total 

100 . 0 % 

57.8%* 

64%* 

94%* 


* All values are based on the total at the maximum volume position. 
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The approximate relative volumes are for a 90° phase difference. For 
the GPU-3 the ratio of maximum volume to minimum volume Is approximately 
1.73. The relative volume and the associated compression and expansion 
piston positions are Illustrated schematically in Figure 5. The sine 
curves represent the position of the pistons with the compression piston 
leading the expansion piston by 90°. The vertical distance between the 
sine curves Is approximately proportional to the gas volume. Note the 
change In volume with crank rotation as Illustrated by the cyclical change 
as the time varies along the horizontal axis. The relative magnitudes of 
significant volumes are Indicated on the figure. 

The curves of Figure 4 Indicate typical gas temperature variation 
(AT/At) due to the compression and expansion effects. The temperature 

C“6 

of the gas is affected by factors other than compression and expansion during 
the cycle. Such factors are not Included here. (See Figures 7, 10 and 11 
and the related discussion.) The curves of Figure 4 represent the solution 
of the isentroplc compression expansion equation 

(AT g /At) ce i (1-y) (T 1 /V 1 )A p R c a)(sin(ti)t+^)-sinu)t) (38) 

The skewed shape of the curves Is primarily due to the dependence on the gas 
temperature which has a skewed Irregular shape as shown In Figure 7. The 
shape of the gas temperature curves is the result of the complex phenomena 
including convective heat transfer , periodic motion through the heater* 
regenerator and cooler regions* and compression and expansion. 

The curves of Figure 6 illustrate the typical cyclic gas pressure 
characteristic in the regenerator. Note that the pressure curves have a 
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Piston Position — Relative Volume AT/At) , °K/sec 
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Schematic Diagram of Relative Volume and Position of Pistons 
(Sine curves represent piston positions , distance between 
curves is approximately proportional to typical gas volume) 
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Figure 7. 


Typical Gas Temperature 


vs Crank Position (time) Curves 


regular shape but are distorted from sine waves with broad "valleys" and 
smooth but narrower peak regions. 

The cited curves are referred to as typical because they illustrate 
general characteristics. With minor modifications and scale changes they 
could be used for computations with different input data. 

Figures 8a and 9a are pressure-volume plots for hydrogen and for 
helium as the working gas. The figures include separate closed curves for 
the compression and the expansion cylinders plotted in a conventional manner 
as* for example* in Ash * ^ . The pairs of closed curves are plotted with the 
compression cylinder volume for the narrower curve, and with the expansion 
cylinder volume for the "fatter" curve. This author questions the signifi- 
cance of such curves because in a Stirling cycle engine the cylinders are 
connected without valving by the cooler, regenerator and heater. The curves 
do not represent separate portions of a cycle as do the p-v diagrams for an 
internal combustion engine. 

The p-v curves of Figures 8b and 9b are believed to be more significant. 
They are plotted as pressure versus volume curves with the volume varying 
from the minimum volume (the volume not displaced plus the minimum volume of 
gas in one or the other or both cylinders) to the maximum volume. Figure 5, 
previously discussed, illustrates these points. As discussed in a previous 
section of this report, the maximum volume for the assumed configuration 
occurs at 135° for the compression piston, and the minimum volume occurs at 
the 315° compression piston position. For the plots of Figures 8b and 9b 
the abscissa is labeled in terms of the decimal fraction of maximum volume. 
For this type of p-v diagram, during the expansion portion from minimum to 
maximum clockwise (the upper curve) the lower of the compression cylinder 
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Pressure N/m x 10 N/m' 


Cylinder Gas Volume/Cyl Displacement 
a. Separate Comp Cyl & Exp Cyl Plots 
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Gas Volume/Max Gas Volume 
b. Plot from Kin Vol to Max Vol 


Figure 8. Typical Pressure-Volume Diagrams for Stirling Cycle 
with Helium at 30 x 1C N/m^ Mean Pressure 


Cylinder Gas Volume/Cyl Displacement 
a. Separate Comp Cyl & Exp Cyl Plots 


Gas Volume/Max Gas Volume 
b. Plot from Min Vol to Max Vol 


Figure 9. Typical Pressure-Volume Diagrams for Stirling Cycle 

5 2 

with Hydrogen at 45 x 10 N/m Mean Pressure 





pressure and the expansion cylinder pressure should be used. During the 
compression portion (the lover part of the closed curve) the higher of the 
compression cylinder pressure and the expansion cylinder pressure should be 
plotted. Using the above logic the specific pressure point plotted will 
depend on the direction of the gas flow. The pressure-volume area (Indicated 
work) will then exclude pressure drop due to fluid flow resistance. 

A typical gas temperature versus position plot is presented in Figure 10. 
This data is for conditions selected to be similar to those used by NASA 
for experimental runs with the results listed in reference**. The plot drawn 
by DEEPP is from computer output data as stored in CEDAT 198. 

Figure 11 is a "softened" gas temperature versus position plot drawn to 
provide a clearer picture of the temperature pattern as it varies in the flow 
region with position and with time. The system has been softened by reducing 
the effective heat transfer in the regenerator to about one-eighth of the 
computed typical value. This computer plot was made from data generated by 
the computer and stored in file CEDAT 63. 

Power versus RPM curves are presented in Figure 12a. These curves for 
the engine operating with hydrogen have been plotted from data in CEDAT 194 
to CEDAT 199. A plot of NASA experimental data as reported in reference"* 
is Included in Figure 12a. 
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FIGURE 11 GAS TEMPERATURE VS. POSITION 

(Softened Syetea) 
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The curves of Figure 12b are typical plots of power versus RPM 
for the assumed Stirling cycle engine configuration using helium as 
the working gas* In Figures 12a and 12b the dasheu portion of the curves 
are based on runs with some variation in conditions. The dashed exten- 
sions while not simple extrapolations are not plotted directly from data 
and operating conditions fully consistent with the sold line portions of 
the curves. 



a. Hydrogen 



Figure 12. Power vs Engine Speed for Various Mean Pressures 
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Effects of Changed Operating Conditions 


9 

i 

The most significant factor affecting computed power output is the 
gas used. Ht^ium and hydrogen were considered. The computations indicated 
that for similar conditions the maximum power output is at considerably 
higher rotational speeds when operating with hydrogen. This is primarily 
due to the lower viscosity and the resultant lower flow pressure drop. 

Other properties of hydrogen also have favorable effects. The power output 
at any given speed la higher for hydrogen than for helium with similar 
assumed conditions. 

Increasing the gas pressure has a direct effect on the computed power 
output. At relatively low pressures the effect appears to be greater than 
a simple linear effect; however the computations did not take into account 
leakage or other problems more pronounced with higher pressures. The effect 
of cooler and heater temperatures is as expected; however only a few variations 
of these parameters were made. Detailed parameter variation studies have not 
been made for other factors. One run was made with much lower temperatures 
and changing from expansion piston phase lag to expansion piston phase leaa. 

This simulated a Stirling cycle heat pump application. The results were not 
analyzed in detail but they were consistent with expectations. 

Suggested Extensions and Modifications of this Research 

It is interesting to note that the temperature variations as computed 
by CDEE for the cooler or compression region are somewhat less than they 

e 

would be for isectroplc compression and expansion while the variations are 
somewhat more than the lsentroplc variations in the heater or expansion 
region. Typical values of the variations are as follows; 
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1. With hydrogen as the working g?s 

a. Cooler AT (CDEE computation) ^ 35* C 

s 

AT (Isentropic) • 45* C 

8 

b. Heater AT (CDEE computation) ± 180* C 

8 

AT (Isentropic) ^ 150* C 

8 

2. With helium as the working gas 

a. Cooler AT (CDEE computation) • 50* C 

8 

AT (Isentropic) ^ 65* C 

8 

b. Heater AT (CDEE computation • 250* C 

8 

ATg (Isentropic) 210* C 

The typical isentropic value' cited are primarily dependent on the comprefsion 
ratio and the value of y ■ C^/Cy, as the cooler and heater temperature levels 
will not vary drastically. Typically for this study T - 15* C and T^ » 704* C. 
The values computed by CDEE are dependent on a number of parameters; the most 
significant ones are compression ratio, temperature levels, heat transfer area, 
the value of y, and other gas properties values. 

The variation from the isentropic case can be accounted for by the heat 
transfer effects. A check of the computed values Indicates that effective 
polytropic "n” values with hydrogen as the working gas are approximately 
1.3 for the cooler and 1.52 for the heater. These values can be compared with 
a ‘V value of 1.4. When helium is used as the working gas the poly tropic "n" 
values are approximately 1.48 for the cooler and 1.85 for the heater, compared 
with a "y" value of 1.67. The regenerator between the cooler and heater is 
relatively close to an Isothermal region with rather small temperature changes 
with time. There is, of course, a large temperature gradient (variation 
with position) in the regenerator. 
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It night be possible to develop $ simple, effective and reasonably 
accurate Stirling cycle analysis technique based on a generalisation and 
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extension of the observations In the preceding paragraph. This would probably 
entail the deternlnatlon of effective polytropic exponents based on design 
factors and correlation with experimental data, followed by thermodynamic , 
fluid mechanics and general mechanical analysis techniques. 

A modification In CDEE that might be worthwhile was tried but not carried 
out to a conclusive evaluation. The intent was to provide for a coarse 
initial coaputatlon by bypassing the computer adjusted time increaent routine. 
This routine aay impose higher accuracy criteria than is justified for scae 
parts of the computation. The intent was to provide an optional QUICK logic 
routine. Perhaps a simple relaxation of intermediate accuracy requirements 
would be a more fruitful approach. This could be accomplished by a simple 
logic addition to permit variable accuracy requirements with a lower accuracy 
for an initial coaputatlon. 

The program would be suitable for a detailed parametric study, and with 
some refinements for system design and optimisation studies. 

A study of transient effects and system control studies could be carried 
out with program simplification so that particularly pertinent data would be 
obtained with minimum computation. The computer program could be modified to 
provide for an iteration on the computed gas density, temperature and pressure 
values. Such a modification could permit larger time and position increments 
and provide greater accuracy. 


- 38 - 


NOMENCLATURE 


t 


Aj area of flow 

A 1 convective heat transfer area 

C c specific heat at constant pressure 

C p specific heat at constant volume 

D v effective diameter 

f gas flow friction factor 

F viscous force 

G V acceleration of gravity 

h convective heat transfer coefficient 

Q 

I increment number 

K pressure coefficient 

K p flow velocity factor 

L effective length 

LC cooler length 

LH heater length 

LR regenerator length 

NINCR number of regenerator increments 
P pressure 

R increment length ratio 

T temperature (°K) 

t time 

U velocity 

U velocity of compression piston 

u velocity of expansion piston 

V volume 

W work 

x position 

x distance to compression piston 

x g distance to expansion piston 

a alpha characteristic line (or point) 

8 beta characteristic line (or point) 

y ratio of specif heats (C /C ) 

A small increment p v 

p density 

p viscosity 

3 partial derivative __ 

7 differential operator (i3/3x+j 3/3y+ka/az) 

w angular velocity (radian /sec) 


Subscripts 


c compression end 

e expansion end 

g gas 

h heater 

i ideal 

m mesh 

r regenerator 

t partial differentiation with respect to t 

x partial differentiation with respect to x 
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